home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / cramer.pro < prev    next >
Text File  |  1997-07-08  |  3KB  |  106 lines
  1. ;$Id: cramer.pro,v 1.6 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       CRAMER
  8. ;
  9. ; PURPOSE:
  10. ;       This function solves an n by n linear system of equations 
  11. ;       using Cramer's rule.
  12. ;
  13. ; CATEGORY:
  14. ;       Linear Algebra.
  15. ;
  16. ; CALLING SEQUENCE:
  17. ;       Result = CRAMER(A, B)
  18. ;
  19. ; INPUTS:
  20. ;       A:      An N by N array of type: float, or double.
  21. ;
  22. ;       B:      An N-element vector of type: float, or double.
  23. ;
  24. ; KEYWORD PARAMETERS:
  25. ;       DOUBLE: If set to a non-zero value, computations are done in
  26. ;               double precision arithmetic.
  27. ;
  28. ;       ZERO:   Use this keyword to set the value of floating-point
  29. ;               zero. A floating-point zero on the main diagonal of
  30. ;               a triangular matrix results in a zero determinant.
  31. ;               A zero determinant results in a 'Singular matrix'
  32. ;               error and stops the execution of CRAMER.PRO.
  33. ;               For single-precision inputs, the default value is 
  34. ;               1.0e-6. For double-precision inputs, the default value 
  35. ;               is 1.0e-12.
  36. ;
  37. ; EXAMPLE:
  38. ;       Define an array (a).
  39. ;         a = [[ 2.0,  1.0,  1.0], $
  40. ;              [ 4.0, -6.0,  0.0], $
  41. ;              [-2.0,  7.0,  2.0]]
  42. ;
  43. ;       And right-side vector (b).
  44. ;         b = [3.0, 10.0, -5.0]
  45. ;
  46. ;       Compute the solution of the system, ax = b.
  47. ;         result = cramer(a, b)
  48. ;
  49. ; PROCEDURE:
  50. ;       CRAMER.PRO uses ratios of column-wise permutations of the array (a)
  51. ;       to calculate the solution vector (x) of the linear system, ax = b.
  52. ;
  53. ; REFERENCE:
  54. ;       ADVANCED ENGINEERING MATHEMATICS (seventh edition)
  55. ;       Erwin Kreyszig
  56. ;       ISBN 0-471-55380-8
  57. ;
  58. ; MODIFICATION HISTORY:
  59. ;       Written by:  GGS, RSI, February 1994
  60. ;       Modified:    GGS, RSI, November 1994
  61. ;                    Added support for double precision results.
  62. ;       Modified:    GGS, RSI, April 1996
  63. ;                    Modified keyword checking and use of double precision.
  64. ;-
  65.  
  66. FUNCTION Cramer, A, B, Double = Double, Zero = Zero
  67.  
  68.   ON_ERROR, 2  ;Return to caller if error occurs.
  69.  
  70.   TypeA = SIZE(A)
  71.   TypeB = SIZE(B)
  72.  
  73.   if TypeA[1] ne TypeA[2] then $
  74.     MESSAGE, "Input array must be square."
  75.  
  76.   if TypeA[3] ne 4 and TypeA[3] ne 5 then $
  77.     MESSAGE, "Input array must be float or double."
  78.  
  79.   ;If the DOUBLE keyword is not set then the internal precision and
  80.   ;result are identical to the type of input.
  81.   if N_ELEMENTS(Double) eq 0 then $
  82.     Double = (TypeA[TypeA[0]+1] eq 5 or TypeB[TypeB[0]+1] eq 5) 
  83.  
  84.   if N_ELEMENTS(Zero) eq 0 and Double eq 0 then $
  85.     Zero = 1.0e-6  ;Single-precision zero.
  86.   if N_ELEMENTS(Zero) eq 0 and Double ne 0 then $
  87.     Zero = 1.0d-12 ;Double-precision zero.
  88.  
  89.   DetermA = DETERM(A, Double = Double, Zero = Zero, /Check)
  90.   if DetermA eq 0 then MESSAGE, "Input array is singular."
  91.  
  92.   if Double eq 0 then xOut = FLTARR(TypeA[1]) $ 
  93.   else xOut = DBLARR(TypeA[1])
  94.  
  95.   for k = 0, TypeA[1]-1 do begin
  96.     ColumnK = A[k,*] ;Save the Kth column of a.
  97.     a[k,*] = B ;Permute the Kth column of A with B.
  98.                ;Solve for the Kth component of the solution xOut
  99.     xOut[k] = DETERM(A, Double = Double, Zero = Zero) / DetermA
  100.     a[k,*] = ColumnK ;Restore A to its original state.
  101.   endfor
  102.  
  103.   RETURN, xOut
  104.  
  105. END
  106.